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ABSTRACT 

Carrier diffusion and thermal conduction play a fundamental role in the operation of high-power, broad-area semi- 
conductor lasers. Restricted geometry, high pumping level and dynamic instability lead to inhomogeneous spatial 
distribution of plasma density, temperature, as well as light field, due to strong light-matter interaction. Thus, 
modeling and simulation of such optoelectronic devices rely on detailed descriptions of carrier dynamics and energy 
transport in the system. 

A self-consistent description of lasing and heating in large-aperture, inhomogeneous edge- or surface-emitting 
lasers (VCSELs) require coupled diffusion equations for carrier density and temperature. In this paper, we derive 
such equations from the Boltzmann transport equation for the carrier distributions. The derived self- and mutual- 
diffusion coefficients are in general nonlinear functions of carrier density and temperature including many-body 
interactions. We study the effects of many-body interactions on these coefficients, as well as the nonlinearity of 
these coefficients for large-area VCSELs. The effects of mutual diffusions on carrier and temperature distributions 
in gain-guided VCSELs will be also presented. 
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1. INTRODUCTION 

In order to achieve higher semiconductor laser output power, two design rules are implemented on chip level -broad 
active area and arrayed structure. As it turns out, the right optical phase locking is needed to achieve efficient 
generation of stable fax-field laser output . 1 Both changes in the carrier density and temperature modulate the optical 
phase through carrier-induced phase shift, and the effect increases in importance as the active region enlarges. At 
the same time, requirements for more stringent control of the transverse-mode dynamics arise since mode spacing 
decreases and more modes become dynamically active, as a result of this enlargement and increase in current injection 
]pvpI 2 jp filamentation and more complex pattern formation and dynamics occur as laser operates at high 

power level . 3 ~ 5 On the other hand, theoretical challenges faced in describing high-power, broad-area semiconductor 
lasers at device-physics level are to fully understand the interplay of carrier dynamics and laser field distribution as 
the ultimate output quality of the laser beam depends on the near-field pattern and spatial phase relation. Since 
stimulated interaction affects the former and carrier dynamics influences the latter, plus phase modulation and 
locking become a sensitive function of the plasma temperature as device size increases, in addition to the necessity 
of a self-consistent framework for the description of an active device on the thermodynamic level, it is thus natural 
to take the plasma temperature and nonlinearity of the physical system into account. Previous attempts to model 
such a dynamic system involve both empirical method 2,4 and first-principle treatment . 6 The advantage of the former 
approach is easy manageability and computational economy, while the latter features fundamental rigor and generic 
inclusiveness. In contrast, a hydrodynamic model starts from first principle, retains part of its rigor and inclusiveness, 
but allows greater numerical efficiency. The purpose of the present paper is to summarize such a model we have 
developed from Boltzmann transport equations (BTEs), in parallel with the semiconductor Bloch equation (SBE), 
discuss its impact on simulations, and present numerical results for gain-guided vertical-cavity surface-emitting lasers 
(VCSELs). 

Further author information: 

J.Z.L.: E-mail: jianzhngOnas.nasa.gov 
S.H.C.: E-mail: cheungOnas.nasa.gov 
C.Z.N.: E-mail: cnmg@nas.nasa.gov 



This paper is organized as follows. In the next section, theoretical results will be summarized, and expressions 
for self- and mutual-diffusion coefficients for carrier density and plasma temperature will be given. Then, numerical 
results and discussions are presented, focusing on several major effects in our model. Finally, we will sum up and 
make some observations for future work. 


2. THEORY 

Now let us summarize the hydrodynamic model we have developed for the description of carrier dynamics in the 
active region of a semiconductor quantum well (QW) laser. We start from the Boltzmann transport equations 7 
(BTEs) for the non-equilibrium carrier distributions n a (k ? r) for electrons (a = e) and holes (a = h) in parabolic 
energy subbands, where k and r are all two-dimensional vectors. Only one subband each is treated for electrons and 
holes. Collisional contributions accounted for in the BTEs are the dominant ones: carricr-LO (longitudinal optical) 
phonon scattering, intraband (both a-a and eh) carrier scattering. As a side note, carrier scattering is explicitly 
retained before we make any assumptions and approximations, which is where our approach differ from others. 6,8 
This unique feature allows a natural and convenient handling of ambipolar transport regime, as it turns out, in the 
development of our model. Moment equation method 9 is used for derivation of the model. The nth-order moment 
and associated current are defined by summing over all degrees of freedom the non-equilibrium distribution function 
with a corresponding weight function F% as below, 

r n {f) = f £ ; /“(f) = | £ q F“nf , (1) 

k k 

where S is the active region area, F* is 1, hk , and h 2 k 2 /2m a , for n = 0,1,2, respectively, and v~ — hk/m a . The first 
three moments and currents are conventionally denoted as: density N a = , momentum P Q = ip f , energy E a = 

density current = Jg, momentum current Jg = , energy current Jg = Jg . Then, it is straightforward to show 

that the first three moment equations can be written as 

8 t N a + d?-J% = B!Z fJ ( 2 ) 

d t P a + d f • r p + N°d? (<5e“ + q a $) = + d t P a U + d t P«\ LO , (3) 

dt E a + df • Jg T df (6e« + q a $) - = R% + d t E% h + d t E a \ LO , (4) 

for electrons (a = e) and holes (a = h). The right-hand-side terms are the results of summing up the corresponding 
generation-recombination (g-r) terms in BTEs. Both intraband (eh) and carrier-LO phonon ( LO ) scatterings couple 
different moment equations. 

It is worth noting that Eqs. (2-4) are exact, but not in closed form — the so-called hierarchy problem. As an 
approximate solution, we assume quasi-equilibrium for the electron-hole plasma, as in fluid dynamics where a drifted 
Maxwell distribution is assumed. 10 We point out that the premise of quasi-equilibrium, upon which further develop- 
ment of our model is hinged, is well established in semiconductor laser theory. 11-13 It means that electrons and holes 
in the plasma, driven out of thermal balance by the laser field and external pumping, are individually described by 
equilibrium distributions of their subsystems in an inertia frame of reference. The physical mechanism behind this 
assumption is the ultrafast carrier intraband (a-a) scattering on the femtosecond time scale. 14 The quasi-equilibrium 
is characterized by the drifted Fermi-Dirac (DFD) distributions, = f a (k — £g) = 1/{1 -b exp[/3 a (e^ - Ma)]} 5 

where kp is the drift wave vector and is the chemical potential. The drift wavevector is related to the first order 

moment and the chemical potential is given by carrier density and temperature. Thus, a total of six parameters for 
the electron and hole distribution functions uniquely relate to the six moment variables. With the aid of the known 
functional form of the DFDs, the right-hand-side terms in Eqs. (2-4) therefore can be expressed in terms of carrier den- 
sities, drift wavevectors, and temperatures. With all this and some symmetric considerations, the hierarchy problem 
for n a (A;,r) is solved — the equations are closed. In contrast, the problem for the interband polarization p(fc, f) needs 
more elaborate treatment 15 which will not be discussed here, as we focus on carrier transport in this paper. Now the 
currents, using the DFD functions, of various orders can be easily shown to depend on the moments after defining the 


thermal energy, by subtracting the kinetic energy from the energy, according to W a — E a — P a • P a /2m a N a , where 
pa _ Therefore, the relations between currents and moments allow us arriving at a closed set of equations 

for {N a ,P Q \W a }, a—c\h . Furthermore, since the thermal energy is a function of carrier density and temperature, 
we have the liberty to use another equivalent set of variables, {N Q ,P a ,7b} , a=e|/i, in which we will derive the final 
form of equations for plasma transport. 

The resulting equations for {N a y P a i W a } arc derived as follows, after expressing the currents, the g-r terms, and 
collisional terms in the moments, 
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where - R% + R%{P a ) 2 /2m a {N a ) 2 -P a P^/m a N a , d,W a \ eh = d t E a \ eh -d t P a \ eh P a /m a N a , and d,W Q \ LO = 

d,E a \ L o-d t P a \ L o * P a /m a N a . Equations (5-7) are the general forms of the moment equations describing carrier 
dynamics in a semiconductor laser, and should be solved together with Poisson equation, the polarization equations, 
and Maxwell equation. Terms {jR^ , d t W a \ € h 7 d t W a \LO } differ by nonlinear terms in P a from {R % , d t E a \ eh , 
d t E a \Lo}- As it turns out, P a ’s are proportional to the space gradients, so the differences can be neglected in case 
of weak inhomogeneity. Furthermore, the set of equations above can be simplified by eliminating the momentum 
equations (6) by resorting to the adiabatic approximation. 8 First, for weakly inhomogeneous systems, we can ignore 
the nonlinear terms in F a , ^though numerically they pos_e no problem. Then, the scattering terms can be linearized 
in the first moment, as dtP a \ e h — —7 e hm r (P a /m a — P& /mp) and dtP a \LO — ~7. w ^ere a /?, 7 e h is the 
momentum relaxation rate due to intrabarid carrier scattering, m r =7n e 77i/j/(m e + m/i) i$ the reduced mass, and 7 i Q 
is the momentum relaxation rate due to carrier-LO phonon scattering. As a result, the solutions for the momenta 
are given as 


p Q = Eg [d?W a + N a d f (fc Q + g°*)] _ 7 lo^a ± rn 0 ) [ OfW + N a _d? ^ 

+Tn ^ + 7eh(m a 72o + iLO^Loi™'* + m p) + 7efc(™a7i 0 + m 0~<Lo) ’ 

which has two terms with distinct physical meanings — the first one, due to intraband (eh) scattering, equilibrates 
the two carrier subsystems, and the second one equilibrates each carrier subsystem with the LO phonons. 

Now let us consider some limiting case that allows a one-component description of the plasma dynamics, i.e 
the ambipolar regime. We have mentioned that the ultrafast intraband (a-a) scattering in the femtosecond range 
is the physical mechanism that leads to the establishment of quasi-equilibrium within each subband. At the same 
time, intraband (eh) scattering, which is on the same time scale, 16 * 17 requires self-consistent treatment within the 
framework of quasi-equilibrium, or we need to consider dynamic correlation between electrons and holes imposed by 
the intraband (eh) scattering. Under the premise of quasi-equilibrium, detailed balance (DB) requires 15 that T e = Th, 
and hkp/m e — hkp/rrih, which is intuitively clear since hkp/m a = Up is the drift velocity of the subsystem. These 
two conditions are equivalent to the ambipolar approximation, since they will reduce the original two-component 
problem to a one-component one. However, to derive the ambipolar diffusion coefficients, and understand when and 
why the ambipolar approximation is valid, we need to look at this issue in more details. Specifically, we consider 
the dynamics around the DB state by looking at the equations for momenta and energies with scattering terms 
linearized around the DB state. As we see from Eq. (8), the DB condition for the drift velocity Up is, in general, not 
valid. Nevertheless, when q < 7^, the second term in Eq. (8) can be neglected, and it leads to the conclusion of 
u e D = Up = #, if N e = N h — N. Thus, we have, in the ambipolar regime, 

7 ehL(P /a + W ^ 

N 7go'Yfo( m « + m p) + 'Tehimalio + m 0^Lo)\ 


( 9 ) 



Therefore, validity of the ambipolar diffusion approximation (ADA) for two components of unequal masses demands 
that the intraband (inter-species) scattering be much faster than scattering between LO phonons and each component 
of the plasma (electrons and holes). Since LO phonon scattering depends sensitively on temperature and the intraband 
scattering on density, it is clear that ADA will no longer be valid at high plasma temperature and relatively low carrier 
density. For now, let us continue the discussion of the limiting case of the ambipolar regime. The hydrodynamic 
equality of drift velocities of the carrier subsystems means that N e — N h — N will be conserved if we prepare 
the system neutrally initially, as indicated from the equation of continuity, Eq. (2). Let this be true, then, charge 
neutrality eliminates the need for Poisson equation, and reduces the density equations for electrons and holes into 
one for the carrier density N, that is 


dfN + d r ' • J N = Rjv . (10) 

In the same spirit, we can show from Eq. (4) that T e = 7/ t = T p , as was found in Ref. 18 in case of a weak electric 
field. This is exactly what the DB condition in quasi-equilibrium requires. We comment that the equalities between 
the two temperatures and the two drift momenta are the consequences of the quasi-equilibrium assumption when 
the phonon scattering is much weaker than the intraband scattering. Using the total plasma energy in place of the 
second moment and Eq. (8), the equation for the second moment is given as follows, 

d t W 4- d? ■ Jw ~ Rw + O t W\Lo > (11) 

where W = W e 4- W h , Jw = 2 Wu, Re — R% + R%, d t W\io = d t E e \LO 4- dtE h \LO- In the above equation, weak 
inhomogeneity has been assumed. To convert the above equation into plasma temperature, the functional dependence 
of the thermal energy on carrier density and plasma temperature, W = W{N,T P ) is used, and the equation is found 
as, 


d t T p + d?- Jt + Jn ' Jn ~ dfjw * Jw — jw ( Rw 4- d t W \ LO ) — JnRn ? (12) 

v 

where Jr — jwJw — JnJn> jw = 1 fdjW, jyv = jw d^W . 

The derived equations above, Eqs. (10,12), include many-body corrections and are applicable to neutral electron- 
hole plasma in a semiconductor quantum well. Using Eq. (9), the diffusion coefficients can be determined in terms of 
the macroscopic variables and microscopic quantities for the plasma. First, the currents for density and temperature 
are given as Jyv = Nu and Jt = (2 jwW/N — Jjv)//v» where Jw = 2 Wu is used. Then, in terms of the gradients of 
plasma density and temperature, the coefficients are defined according to 


Jn = —DNNdfN — DNrdfTp , (13) 

J T = —DrNdfN — DjjdfTp . (14) 

Finally, after expressing the drift velocity, Eq. (9), in terms of the gradients of density and temperature, the diffusion 
coefficients are straightforwardly obtained as 
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where W is the total thermal energy of the plasma, Se g = Se e 4- 5e h is the total Fock exchange correction to the 
bandgap, and T u = {ll 0 lLo( m e + m h) + lehi^eJlo + m hJLo)]/jeh* Before moving on, we introduce the following 
terminology for the diffusion coefficients for future discussions. The coefficient is dubbed self-diffusion when it 
connects the gradient of a variable to its current, and mutual- diffusion otherwise. In addition, variable name is used 
to label the coefficient that relates to the variable gradient. For example, Dtn is called mutual- diffusion density 
coefficient accordingly. 


To complete this section, we write down the macroscopic polarization and laser field equations, which arc treated 
as in Ref. 19. The polarization is decomposed to two parts, electronic and background, as P = Po + Pi* where the 
background contribution is given as Po — to^oXo (N, 7j)£, and the electronic part is dynamically described by 

O t Pi = , ( 17 ) 

where the parameters Xo, the effective background susceptibility, Ti, the gain bandwidth, uf i, and A\ y the Lorentzian 
oscillator strength, are fitted to microscopically computed values as a function of carrier density N , plasma temper- 
ature r p , and lattice temperature T \ , 19 The laser field is described, after integrating over the assumed longitudinal 
mode distribution, by 


d t £ = 


+ 

2K r 


iVgKT 

2€ 0 e b 


P-k£ . 


( 18 ) 


Standard notations are used in the above equations, for example, u c is the central frequency, and T is the optical 
mode confinement factor. Finally, the generation-recombination terms in Eqs. (10) and (12) are given as below, 


r n = -7 N n + _ ^S(P'£) , 

e 4 n 

R w = E g - lV g 9 [(hAu c - ihj p )P m £ - » 7A ,ft0 t P*£] , 

where A uj c = u c ~ is the detuning, J is the injection current with a pumping profile and an quantum efficiency 
77, L m is the aggregate active region width, A E g is the bandgap offset, and N s is scaling constant of 10 12 cm" 2 . In 
the equation for Rw, we neglect the generation-recombination contributions and the last term is derived in Ref. 20. 
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Figure 1. Fundamental mode results with (open-circled solid line) and without (x-symbol dashed line) plasma 
heating effects for a gain-guided circular InGaAs VCSEL. Shown here arc the cross-sectional distributions through 
the center. 


3. SIMULATION RESULTS AND DISCUSSIONS 

The simulations are run for a gain-guided InGaAs/GaAs VCSEL operating at u) c of 980 nm. Injection current profile 
is circular with a diameter of 7.5 fim. The cavity length L is 144 nm and L m is 36 nm, which yields a T of 0.25. 
The other VCSEL parameters and simulation details can be found in Ref. 21. Except the carrier lifetime (1/7jv) of 
2.5 ns, other rates are computed microscopically and fitted as a function of N> T p , and TJ. Furthermore, we use an 



injection current J of 1.4D6 KA/cm 2 , which is about 18% above the threshold current. And the lattice temperature 
Ti is assumed to be a constant of 295 K throughout our simulations. 

First, we demonstrate the necessity of including plasma temperature in order to accurately predict lasing per- 
formance of the VCSEL. As shown in Fig. 1, plasma heating (middle panel) leads to roughly 40% reduction in the 
ncar-ficld laser intensity (right panel) and a weaker spatial hole burning effect (left panel). The shoulder structure in 
the plasma temperature distribution is attributed to the restrictive current pumping profile which is disk-like with 
a diameter of 7.5 fim constant pumping region in the middle and a smooth continuation area of total diameter of 
13.5 ^m. Beyond this region the pumping current is zero. As a result, inside the middle region, injection heating 
contributes to a higher plasma temperature, while in the continuation area a dip structure develops in the shoulder 
and explanation is given in Subsection 3.3. Qualitatively, plasma heating effects cause gain degradation, carrier 
density modification, weaker electronic susceptibility, and thus smaller laser intensity. 



Figure 2. Numerical results of diffusion coefficients in log-log scale. Insets (in linear XY scale) show the percentage 
corrections to the coefficients due to the many-body exchange interactions. Indicated are the nonlinear dependence 
of the coefficients on carrier density and many-body effects on them. Results at three temperatures are shown: solid 
line — 200 K, dotted line 300 K, dot-dashed line— 400 K. 

3.1. Many-body effects on diffusion coefficients 

We now discuss the many-body effects on the diffusion coefficients, owing to the attractive Coulomb interactions 
between electrons and holes. Exhibited in the insets of Fig. 2, are the percentage change in the coefficients after the 
Fock exchange term in Eqs. (15-16) is taken into account. The change is around 10% at typical III-V lasing density. 
The dependence of the effects on density falls into two groups: Dnn and Dtn vs D^t and Dtt • The former shows 
a maximum below the density of 10 12 cm -2 and weak dependence at high density, while the latter i ncrease s with 
density. The difference derives from the dependence of the bandgap renormalization (BGR) term 5z g on density and 
temperature. Our numerical results 22 for quasi-2D plasma in 8 nm CaAs/AlGaAs revealed a weaker than linear 
relation power-law dependence for the BGR term on density in the free carrier regime. Similar results is expected 
for the present material. In addition, at low density side, the BGR term increases with density faster than linear 
dependence since the interaction increases as inter-particle separation is reduced. The overall behavior of the former 
group thus can be understood from these understandings. In comparison, the latter group is easier to understand, 
since in the free carrier regime the magnitude of the BGR term increases with plasma temperature due to less effective 
screening at higher temperature. As a result, the ex change correction increases in magnitude sub-linearly according 
to our numerical results. 22 Furthermore, data at three temperatures are shown: 200 K, 300 K, and 400 K. The figure 
reveals that temperature reduces the many-body effects as well as improves the diffusivity of the plasma, which is 
intuitively comprehensible. As temperature increases, the kinetic energy of a typical particle increases accordingly, 
which leaves the barely changed interaction term less significant and makes the particle easier to diffuse around. This 




is why the k B T factor appears in the famous Einstein relation that connects diffusion coefficient and mobility of the 
particle. 
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Figure 3. Nonlinear diffusion effects of diffusion coefficient Djvyv in the InGaAs VCSEL. Shown are temporal 
evolution towards steady state (1 ns) under constant pumping and no laser field and plasma heating. See text for 
curve style usage. 


3.2. Nonlinearity effect 

As demonstrated in Fig. 2, the diffusion coefficients are nonlinear functions of carrier density and plasma temperature. 
In the case of a gain-guided device, the density varies drastically from the active region to the non-lasing area, over 
a few orders of magnitude. Therefore, the nonlinearity in the coefficients is expected to influence the density and 
temperature distributions. For illustrative purpose and simplicity, we only consider the nonlinearity in the coefficient 
Dux and neglect the temperature equation totally. Further, we consider two scenarios under constant current 
injection: without lasing (Fig. 3) and with lasing (Fig. 4). Four cases are studies for comparison in each scenario: (a) 
full nonlinearity— solid line; and three linear cases (b) D^n = 108.24 cm 2 /s — dotted line; (c) 54.12 cm 2 /s dashed 
line; and (d) 23.84 cm 2 /s— dot-dashed line. The values of the coefficient are chosen in the linear cases corresponding 
to different densities and the same temperature of300K: (b) 5xl0 12 cm” 2 ; (c) 2.5xl0 12 cm” 2 ; and (d) 7xl0 n cm” 2 . 



Figure 4. Nonlinear effects of diffusion coefficient Djvjv in the InGaAs VCSEL. Same as Fig. 3 except that laser 
field is included. Note that, in the linear case with the largest coefficient (case b: dotted line), no lasing is predicted. 






We discuss Fig. 3, the scenario with only carrier density diffusion, first. Shown are four snapshots (120 ps, 300 
ps, GOO ps, and 1 ns) of the transient dynamics of the VCSEL under the same initial condition and evolving towards 
a steady state (1 ns). Note that the injection current is increased to 10.8 KA/cm 2 accidentally here and have no 
influence on our comparison. Clearly, the larger the coefficient, the more spreaded out the distribution, lower the 
center peak, and faster the steady state approached. Consequently, these effects impact on the lasing performance. 
As we look at Fig. 4 now, attention is called to curve b (dotted line), which has the largest coefficient and does not 
lase. The reason is simple: not enough gain available, as indicated in the 1 ns steady state data. The carrier density 
at the peak position is below transparency as compared to curve c (dashed line). Summarily, the nonlinear effects 
of the diffusion coefficients affect (1) the steady state distribution, (2) the transient dynamics, (3) threshold current, 
and thus (4) laser output power and spatial hole burning effect. 
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Figure 5. Mutual-diffusion effects in a gain-guided InGaAs VCSEL. Solid curve with open circles includes mutual- 
diffusion terms, and the dotted curve with x symbol is without those terms. In fundamental mode operation, “no 
discernable effect in the density and near-field distributions of the laser is shown. The dips in the shoulder structure 
of the plasma temperature are gone when there are no mutual diffusions present. 

3.3- Mutual diffusions between density and temperature 

In this study we aim to examine the importance of the mutual-diffusion terms, first introduced by our model, 
in the simulation. We compare two cases: with (open-circled solid line) and without (x- symbol dotted line) the 
mutual-diffusion terms in Eqs. (13 and 14), as they appear in Eqs. (10) and (12). In fig. 5, the carrier density, 
plasma temperature, and the near-field laser intensity are shown. In fundamental mode operation, under moderate 
current pumping levd as we are using in the present simulations, no discernable effects in the density and near-field 
distributions of the laser are found. The reason is mainly that the coefficient values (refer to Fig. 2) tell that the 
density gradient has a larger effect than the temperature gradient. This is corroborated in the temperature data 
(center panel), as neglect of the mutual-diffusion terms changes its distribution and the dips in the shoulder structure 
are removed. This revelation also explains the cause of the dips. Finally, the extraordinary discontinuities in the 
data outside of the shoulder structure in plasma temperature are due to numerical glitches and under investigation. 

4. SUMMARY 

In this paper, we summarize our recently developed hydrodynamic model for semiconductor lasers. The model 
includes plasma heating effects and is derived from first principles. In the case where carrier scattering predominates 
over carrier-phonon scattering, an junbipolar regime is realized and the plasma can be described like a one-component 
physical system. This is the regime in which semiconductor lasers operate normally. We give the expressions for 
diffusion coefficients in this regime and investigate some major effects in our model for a exemplary, gain-guided, 
single-mode VCSEL with an diameter of 7.5 /im. Nonlinear effects of one of the coefficients, JDjvat, on the performance 
of the device under both non-lasing and lasing conditions are studied. It is found that the influence is appreciable in 



accurately predicting the performance of the device. Also pointed out is the necessity of inclusion of plasma heating 
in the simulation of semiconductor lasers. Finally, we study the effect of mutual-diffusion interaction between carrier 
density and plasma temperature, and the effect is found to be negligible. Future studies will consider index-guided, 
larger area VCSELS and arrays. 
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